///////////////////////////////////////////////////////////////////////////////////////////////////
/// \file configuration.h
///
/// \brief   header file for the Configuration and ConfigDiff classes
///
///////////////////////////////////////////////////////////////////////////////////////////////////

#ifndef PCPS_CONFIGURATION_HEADER
#define PCPS_CONFIGURATION_HEADER

#include<src/pcps.h>
#include<src/input.h>
#include<src/random.h>

namespace pcps {

  template<class S> class Configuration;

  ///////////////////////////////////////////////////////////////////////////////////////////////////
  /// \brief a class to hold information about the difference between two configurations
  ///
  ///////////////////////////////////////////////////////////////////////////////////////////////////
  template<class S> class ConfigDiff {

    private:

      /// \brief  sites in the fock vector that are empty in new but occupied in old
      std::vector<int> _newly_empty;

      /// \brief  sites in the fock vector that are occupied in new but empty in old
      std::vector<int> _newly_full;

      // /// \brief  A vector of characters telling which of the vectors in each u-v pair are changing.
      // ///         Allowed characters are u, v, and n (for niether).
      // std::vector<char> _noise_change_codes;

      // /// \brief  For each u-v pair, the difference given by the new u or v minus the old u or v.
      // std::vector<std::vector<S> > _noise_changes;

    public:

      std::vector<int> & newly_empty() { return _newly_empty; }
      std::vector<int> & newly_full() { return _newly_full; }
      // std::vector<char> & noise_change_codes() { return _noise_change_codes; }
      // std::vector<std::vector<S> > & noise_changes() { return _noise_changes; }

      const std::vector<int> & newly_empty() const { return _newly_empty; }
      const std::vector<int> & newly_full() const { return _newly_full; }
      // const std::vector<char> & noise_change_codes() const { return _noise_change_codes; }
      // const std::vector<std::vector<S> > & noise_changes() const { return _noise_changes; }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief default constructor
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      ConfigDiff()
        : _newly_empty(),
          _newly_full()
          //_noise_change_codes(2, 'n'),
          //_noise_changes(2)
      {}

      ConfigDiff(const pcps::Configuration<S> & new_config, const pcps::Configuration<S> & old_config);

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief Constructs the configuration difference object
      ///
      /// \param[in]      userinp  the user's input options
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      ConfigDiff(const pcps::Input & userinp)
        : _newly_empty(),
          _newly_full()
          //_noise_change_codes(2, 'n'),
          //_noise_changes(2)
      {
        //if (userinp.z_spin() != 0.0)
        //  throw pcps::Exception("config difference constructor currently assumes equal number of alpha and beta");
        //for (int i = 0; i < _noise_changes.size(); i++)
        //  _noise_changes.at(i).assign(userinp.nparticles()/2, pcps::zero<S>());
      }

  };

  ///////////////////////////////////////////////////////////////////////////////////////////////////
  /// \brief  a class to hold a configuration
  ///
  ///////////////////////////////////////////////////////////////////////////////////////////////////
  template<class S> class Configuration {

    private:

      /// \brief  the number of alpha sites
      int _nas;

      /// \brief  the number of beta sites
      int _nbs;

      /// \brief  the number of alpha particles
      int _nap;

      /// \brief  the number of beta particles
      int _nbp;

      /// \brief  the number of alpha holes
      int _nau;

      /// \brief  the number of beta holes
      int _nbu;

      /// \brief  the configuration of the fock vector, in ababab... format
      std::vector<int> _fock_config;

      /// \brief  for each alpha site, the position of that site (if occupied) in the list of occupied sites
      std::vector<int> _a_site_to_occ;

      /// \brief  for each beta site, the position of that site (if occupied) in the list of occupied sites
      std::vector<int> _b_site_to_occ;

      /// \brief  for each alpha site, the position of that site (if unoccupied) in the list of unoccupied sites
      std::vector<int> _a_site_to_unocc;

      /// \brief  for each beta site, the position of that site (if unoccupied) in the list of unoccupied sites
      std::vector<int> _b_site_to_unocc;

      /// \brief  the site indeces of the occupied alpha sites
      std::vector<int> _a_occ_to_site;

      /// \brief  the site indeces of the occupied beta sites
      std::vector<int> _b_occ_to_site;

      /// \brief  the site indeces of the unoccupied alpha sites
      std::vector<int> _a_unocc_to_site;

      /// \brief  the site indeces of the unoccupied beta sites
      std::vector<int> _b_unocc_to_site;

      // /// \brief  the number of noise vectors to use
      // int _n_noise_vecs;

      // /// \brief  the value of the u noise vectors used for removing nodes from determinants, agp, etc.
      // std::vector<std::vector<S> > _noise_vecs_u;

      // /// \brief  the value of the v noise vectors used for removing nodes from determinants, agp, etc.
      // std::vector<std::vector<S> > _noise_vecs_v;

    public:

      std::vector<int> & fock_config() { return _fock_config; }
      //std::vector<std::vector<S> > & noise_vecs_u() { return _noise_vecs_u; }
      //std::vector<std::vector<S> > & noise_vecs_v() { return _noise_vecs_v; }
      int & operator[](const int i) { assert( i >= 0 && i < _fock_config.size() ); return _fock_config[i]; }

      const std::vector<int> & fock_config() const { return _fock_config; }
      //const std::vector<std::vector<S> > & noise_vecs_u() const { return _noise_vecs_u; }
      //const std::vector<std::vector<S> > & noise_vecs_v() const { return _noise_vecs_v; }
      int operator[](const int i) const { assert( i >= 0 && i < _fock_config.size() ); return _fock_config[i]; }

      const std::vector<int> & a_site_to_occ   () const { return _a_site_to_occ;   }
      const std::vector<int> & b_site_to_occ   () const { return _b_site_to_occ;   }
      const std::vector<int> & a_site_to_unocc () const { return _a_site_to_unocc; }
      const std::vector<int> & b_site_to_unocc () const { return _b_site_to_unocc; }
      const std::vector<int> & a_occ_to_site   () const { return _a_occ_to_site;   }
      const std::vector<int> & b_occ_to_site   () const { return _b_occ_to_site;   }
      const std::vector<int> & a_unocc_to_site () const { return _a_unocc_to_site; }
      const std::vector<int> & b_unocc_to_site () const { return _b_unocc_to_site; }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief  return a compound index for a pair of sites, one occupied and the other empty
      ///
      /// \param[in]     occ_site     the index of the occupied site (in alpha and beta indexing)
      /// \param[in]     unocc_site   the index of the unoccupied site (in alpha and beta indexing)
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      int occ_unocc_cmpd(int occ_site, int unocc_site) const {

        const int b = occ_site % 2;

        assert( occ_site % 2 == unocc_site % 2 );

        const std::vector<int> & site_to_occ   = ( b == 0 ? _a_site_to_occ   : _b_site_to_occ   );
        const std::vector<int> & site_to_unocc = ( b == 0 ? _a_site_to_unocc : _b_site_to_unocc );
        const int np = ( b == 0 ? _nap : _nbp );
        const int nu = ( b == 0 ? _nau : _nbu );

        occ_site /= 2;
        unocc_site /= 2;

        assert( site_to_occ[occ_site] >= 0 );
        assert( site_to_occ[occ_site] < np );

        int retval;

        if ( occ_site == unocc_site ) {

          retval = site_to_occ[occ_site];

        } else {

          assert( site_to_unocc[unocc_site] >= 0 );
          assert( site_to_unocc[unocc_site] < nu );

          retval = np + nu * site_to_occ[occ_site] + site_to_unocc[unocc_site];

        }

        assert( retval >= 0 );
        assert( retval < np*(nu+1) );

        return retval;

      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief proposes a move from one configuration to another
      ///
      /// \param[in]     userinp      the user's input options
      /// \param[in]     lfg          a random number generator
      /// \param[in]     config       the configuration to move from
      /// \param[out]    new_config   on exit, the configuration to move to
      /// \param[out]    config_diff  on exit, specifies the difference between the configurations
      ///
      /// \return the proposal density ratio, which is the ratio of the probability of proposing the
      ///         current configuration if we are at the new configuration to the probability
      ///         of proposing the new configuration if we are at the current configuration.
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      static double propose_move(const pcps::Input & userinp,
                                 pcps::LaggedFibonacci<> & lfg,
                                 const pcps::Configuration<S> & config,
                                 pcps::Configuration<S> & new_config,
                                 pcps::ConfigDiff<S> & config_diff) {

        // // create a function object for creating new random noise vectors
        // struct NoiseUpdater {
        //   static void update(const double guiding_noise,
        //                      const std::vector<S> & old_vec,
        //                      std::vector<S> & new_vec,
        //                      std::vector<S> & diff_vec) {
        //     const double x = 2.0 * std::sqrt(guiding_noise);
        //     const S y = 0.5 * ( pcps::unity<S>() + pcps::imaginary_unity<S>() );
        //     for (int j = 0; j < new_vec.size(); j++) {
        //       new_vec.at(j) = x * ( pcps::random_number<S>() - y ); // get new vector
        //       diff_vec.at(j) = new_vec.at(j) - old_vec.at(j); // get difference from old vector
        //     }
        //   }
        // };

        // propose a new fock space configuration by moving one alpha or one beta electron
        const int pos1 = lfg() % userinp.nsites();
        const int pos1_occ = config.fock_config().at(pos1);
        assert(pos1_occ == 0 || pos1_occ == 1);
        const int b = pos1 % 2;
        const int n_a_part = ( b == 0 ? userinp.n_elec_a() : userinp.n_elec_b() );
        const int x = lfg() % ( pos1_occ == 1 ? userinp.nsites()/2 - n_a_part : n_a_part );
        int c = 0;
        int j = 0;
        for ( ; j < userinp.nsites()/2; j++)
          if (config.fock_config().at(2*j+b) != pos1_occ)
            if (c++ == x)
              break;
        const int pos2 = 2*j+b;
        const int pos2_occ = config.fock_config().at(pos2);
        assert(c == x+1);
        assert(pos2_occ == 0 || pos2_occ == 1);
        assert(pos1 != pos2);
        assert(pos1_occ != pos2_occ);
        assert(pos1 % 2 == pos2 % 2);

        // define the difference
        config_diff.newly_empty().assign(1, (pos1_occ == 1 ? pos1 : pos2));
        config_diff.newly_full().assign(1, (pos1_occ == 0 ? pos1 : pos2));

        // get the new configuration
        new_config = config;
        new_config.update(config_diff);
        //std::swap(new_config.fock_config().at(pos1), new_config.fock_config().at(pos2));

        // // propose new noise vectors
        // for (int i = 0; i < new_config.noise_vecs_u().size(); i++) {
        //   config_diff.noise_change_codes().at(i) = ( b == 0 ? 'v' : 'u' );
        //   if (config_diff.noise_change_codes().at(i) == 'u')
        //     NoiseUpdater::update( userinp.guiding_noise(), config.noise_vecs_u().at(i),
        //                           new_config.noise_vecs_u().at(i), config_diff.noise_changes().at(i) );
        //   if (config_diff.noise_change_codes().at(i) == 'v')
        //     NoiseUpdater::update( userinp.guiding_noise(), config.noise_vecs_v().at(i),
        //                           new_config.noise_vecs_v().at(i), config_diff.noise_changes().at(i) );
        // }

        // for now, the proposal density is symmetric, so return one
        return 1.0;

      }

//      ///////////////////////////////////////////////////////////////////////////////////////////////////
//      /// \brief Generates random noise vectors
//      ///
//      /// \param[in]      userinp  the user's input options
//      ///
//      ///////////////////////////////////////////////////////////////////////////////////////////////////
//      void generate_noise(const pcps::Input & userinp) {
//
//        if (userinp.z_spin() != 0.0)
//          throw pcps::Exception("noise generator currently assumes equal number of alpha and beta");
//
//        const S y = 0.5 * ( pcps::unity<S>() + pcps::imaginary_unity<S>() );
//        const double x = 2.0 * std::sqrt(userinp.guiding_noise());
//
//        for (int i = 0; i < _noise_vecs_u.size(); i++) {
//          _noise_vecs_u.at(i).resize(userinp.nparticles()/2);
//          _noise_vecs_v.at(i).resize(userinp.nparticles()/2);
//          for (int j = 0; j < _noise_vecs_u.at(i).size(); j++) {
//            _noise_vecs_u.at(i).at(j) = x * ( pcps::random_number<S>() - y );
//            _noise_vecs_v.at(i).at(j) = x * ( pcps::random_number<S>() - y );
//          }
//        }
//
//      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief Constructs the configuration object
      ///
      /// \param[in]      userinp  the user's input options
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      Configuration(const pcps::Input & userinp)
        : _nas(userinp.nsites()/2),
          _nbs(userinp.nsites()/2),
          _nap( pcps::round( userinp.nparticles()/2.0 + userinp.z_spin() ) ),
          _nbp( pcps::round( userinp.nparticles()/2.0 - userinp.z_spin() ) ),
          _nau( _nas - _nap ),
          _nbu( _nbs - _nbp ),
          _fock_config(userinp.nsites()),
          _a_site_to_occ(_nas, -1),
          _b_site_to_occ(_nbs, -1),
          _a_site_to_unocc(_nas, -1),
          _b_site_to_unocc(_nbs, -1),
          _a_occ_to_site(_nap, -1),
          _b_occ_to_site(_nbp, -1),
          _a_unocc_to_site(_nau, -1),
          _b_unocc_to_site(_nbu, -1)
          //_n_noise_vecs(2),
          //_noise_vecs_u(_n_noise_vecs),
          //_noise_vecs_v(_n_noise_vecs)
      {

        if (userinp.z_spin() != 0.0)
          throw pcps::Exception("configuration constructor currently assumes equal number of alpha and beta");

        // set some constants
        const int nsites = userinp.nsites();
        std::vector<int> ab_n_particles(2, 0);
        ab_n_particles.at(0) = userinp.n_elec_a();
        ab_n_particles.at(1) = userinp.n_elec_b();

        // // initialize the random noise vectors
        // this->generate_noise(userinp);

        // if supplied, start in the user's chosen fock space configuration
        if (userinp.init_config().size() == nsites) {

          _fock_config = userinp.init_config();

        // otherwise, generate a random fock space configuration
        } else {

          // remove all particles from the configuration
          _fock_config.assign(nsites, 0);

          // place alpha (b == 0) and beta (b == 1) particles
          for (int b = 0; b <= 1; b++) {
            for (int i = 0; i < ab_n_particles.at(b); i++) {
              int x = pcps::global_lcg() % (nsites/2 - i);
              int c = 0;
              int j = 0;
              for ( ; j < nsites/2; j++)
                if (_fock_config.at(2*j+b) == 0)
                  if (c++ == x)
                    break;
              assert(c == x+1);
              assert(_fock_config.at(2*j+b) == 0);
              _fock_config.at(2*j+b) = 1;
            }
          }

        }

        // initialize occupation/unoccupation/site mappings
        this->init_occ_unocc_site();

      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief  Initializes the occ <--> site and unocc <--> site mappings for the stored configuration
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      void init_occ_unocc_site() {

        for (int b = 0; b < 2; b++) {

          std::vector<int> & site_to_occ   = ( b == 0 ? _a_site_to_occ   : _b_site_to_occ   );
          std::vector<int> & site_to_unocc = ( b == 0 ? _a_site_to_unocc : _b_site_to_unocc );
          std::vector<int> & occ_to_site   = ( b == 0 ? _a_occ_to_site   : _b_occ_to_site   );
          std::vector<int> & unocc_to_site = ( b == 0 ? _a_unocc_to_site : _b_unocc_to_site );
          const int ns = ( b == 0 ? _nas : _nbs );
          const int np = ( b == 0 ? _nap : _nbp );
          const int nu = ( b == 0 ? _nau : _nbu );

          assert( _fock_config.size() == 2*ns );

          int pc = 0; // particle counter
          int hc = 0; // hole counter

          int i = b;
          for ( ; i < _fock_config.size(); i+=2) {

            if ( _fock_config[i] == 1 ) {

              site_to_unocc[i/2] = -1;
              site_to_occ[i/2] = pc;
              occ_to_site[pc] = i/2;
              pc++;

            } else {

              site_to_occ[i/2] = -1;
              site_to_unocc[i/2] = hc;
              unocc_to_site[hc] = i/2;
              hc++;

            }

          }

          assert( i/2 == ns );
          assert(  pc == np );
          assert(  hc == nu );

        }

      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief  Updates the object to incorperate the supplied difference.
      ///
      /// \param[in]     diff  specifies the difference between the configuration to update to
      ///                      and the current configuration
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      void update(const pcps::ConfigDiff<S> & diff) {

        assert( diff.newly_full().size() == diff.newly_empty().size() );

        for (int i = 0; i < diff.newly_full().size(); i++) {

          int new_occ = diff.newly_full().at(i);
          int old_occ = diff.newly_empty().at(i);

          assert( new_occ % 2 == old_occ % 2 );

          assert( _fock_config.at(new_occ) == 0 );
          assert( _fock_config.at(old_occ) == 1 );

          _fock_config.at(new_occ) = 1;
          _fock_config.at(old_occ) = 0;

          const int b = new_occ % 2;

          std::vector<int> & site_to_occ   = ( b == 0 ? _a_site_to_occ   : _b_site_to_occ   );
          std::vector<int> & site_to_unocc = ( b == 0 ? _a_site_to_unocc : _b_site_to_unocc );
          std::vector<int> & occ_to_site   = ( b == 0 ? _a_occ_to_site   : _b_occ_to_site   );
          std::vector<int> & unocc_to_site = ( b == 0 ? _a_unocc_to_site : _b_unocc_to_site );
          const int ns = ( b == 0 ? _nas : _nbs );
          const int np = ( b == 0 ? _nap : _nbp );
          const int nu = ( b == 0 ? _nau : _nbu );

          new_occ /= 2;
          old_occ /= 2;

          assert( site_to_occ[old_occ] >= 0 );
          assert( site_to_occ[old_occ] < np );
          assert( occ_to_site[ site_to_occ[old_occ] ] == old_occ );
          occ_to_site[ site_to_occ[old_occ] ] = new_occ;
          site_to_occ[new_occ] = site_to_occ[old_occ];
          site_to_occ[old_occ] = -1;

          assert( site_to_unocc[new_occ] >= 0 );
          assert( site_to_unocc[new_occ] < nu );
          assert( unocc_to_site[ site_to_unocc[new_occ] ] == new_occ );
          unocc_to_site[ site_to_unocc[new_occ] ] = old_occ;
          site_to_unocc[old_occ] = site_to_unocc[new_occ];
          site_to_unocc[new_occ] = -1;

        }

      }

  };

}

#endif
